#此代码中的pattern是针对DC样本统计的，分得很细，共有8种pattern。而双发和健康中只有5种pattern，这部分统计的代码在20201105中。
#二次修改时间是11.15，为了把标准换成FDR<0.1，得出来的结果是一样的，单发样本中pattern一致的仍是22个位点。
setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

i=1
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file$FDR1=p.adjust(file$normal_bayes_pvalue,method = "BH")
file$FDR2=p.adjust(file$tumor_bayes_pvalue,method = "BH")
file=file[file$FDR1<0.1|file$FDR2<0.1,]
file$normal_group="nosig"
file[file$normal_bayes_beta0>0&file$FDR1<0.1,]$normal_group="up"
file[file$normal_bayes_beta0<0&file$FDR1<0.1,]$normal_group="down"
file$tumor_group="nosig"
file[file$tumor_bayes_beta0>0&file$FDR2<0.1,]$tumor_group="up"
file[file$tumor_bayes_beta0<0&file$FDR2<0.1,]$tumor_group="down"
file$pattern=paste0("normal_",file$normal_group,"-tumor_",file$tumor_group)

normal_up_tumor_nosig=file[file$pattern=="normal_up-tumor_nosig",]$unitID
normal_up_tumor_up=file[file$pattern=="normal_up-tumor_up",]$unitID
normal_down_tumor_nosig=file[file$pattern=="normal_down-tumor_nosig",]$unitID
normal_nosig_tumor_up=file[file$pattern=="normal_nosig-tumor_up",]$unitID
normal_down_tumor_up=file[file$pattern=="normal_down-tumor_up",]$unitID
normal_up_tumor_down=file[file$pattern=="normal_up-tumor_down",]$unitID
normal_down_tumor_down=file[file$pattern=="normal_down-tumor_down",]$unitID
normal_nosig_tumor_down=file[file$pattern=="normal_nosig-tumor_down",]$unitID
###单发样本取交集，upup有14个交集，downdown8个，其余的没有交集。
for (i in 2:6) {
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file$FDR1=p.adjust(file$normal_bayes_pvalue,method = "BH")
file$FDR2=p.adjust(file$tumor_bayes_pvalue,method = "BH")
file=file[file$FDR1<0.1|file$FDR2<0.1,]
file$normal_group="nosig"
file[file$normal_bayes_beta0>0&file$FDR1<0.1,]$normal_group="up"
file[file$normal_bayes_beta0<0&file$FDR1<0.1,]$normal_group="down"
file$tumor_group="nosig"
file[file$tumor_bayes_beta0>0&file$FDR2<0.1,]$tumor_group="up"
file[file$tumor_bayes_beta0<0&file$FDR2<0.1,]$tumor_group="down"
file$pattern=paste0("normal_",file$normal_group,"-tumor_",file$tumor_group)
normal_up_tumor_nosig=intersect(normal_up_tumor_nosig,file[file$pattern=="normal_up-tumor_nosig",]$unitID)
  normal_up_tumor_up=intersect(normal_up_tumor_up,file[file$pattern=="normal_up-tumor_up",]$unitID)
  normal_down_tumor_nosig=intersect(normal_down_tumor_nosig,file[file$pattern=="normal_down-tumor_nosig",]$unitID)
  normal_nosig_tumor_up=intersect(normal_nosig_tumor_up,file[file$pattern=="normal_nosig-tumor_up",]$unitID)
  normal_down_tumor_up=intersect(normal_down_tumor_up,file[file$pattern=="normal_down-tumor_up",]$unitID)
  normal_up_tumor_down=intersect(normal_up_tumor_down,file[file$pattern=="normal_up-tumor_down",]$unitID)
  normal_down_tumor_down=intersect(normal_down_tumor_down,file[file$pattern=="normal_down-tumor_down",]$unitID)
  normal_nosig_tumor_down=intersect(normal_nosig_tumor_down,file[file$pattern=="normal_nosig-tumor_down",]$unitID)
}
结果在20201104下的intersect后缀文件中。
write.table(normal_up_tumor_up,"../20201104/normal_up_tumor_up.intersect.txt",quote = F,row.names = F)
write.table(normal_down_tumor_down,"../20201104/normal_down_tumor_down.intersect.txt",quote = F,row.names = F)
normal_up_tumor_up=data.frame(unitID=normal_up_tumor_up,group=rep("normal_up_tumor_up",length(normal_up_tumor_up)))
 normal_down_tumor_down=data.frame(unitID=normal_down_tumor_down,group=rep("normal_down_tumor_down",length(normal_down_tumor_down)))
 rt=rbind(normal_down_tumor_down,normal_up_tumor_up)
 rt1=tidyr::separate(rt,unitID,into=c("chr","pos","ref","alt"),sep=":")
 rt1=data.frame(rt1[,1:2],rt1[2:5])
 write.table(rt1,"intersect.DC.txt",quote=F,row.names = F,sep="\t")
 
 